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A FORTRAN PROGRAM FOR SOLVING SYSTEMS OF COUPLED 
SECOND-ORDER DIFFERENTIAL EQUATIONS WITH TWO-POINT 
BOUNDARY CONDITIONS 


INTRODUCTION 


Systems of second-order ordinary differential equations with two-point 
boundary conditions are invariably encountered in problems of applied physics. 
Lindzen and Kuo (1969) have published an algorithm for solving such systems 
by the method of Gaussian Elimination. The method was found to be of particular 
value in the solution of various thermospheric tidal problems currently under 
investigation at Boston College and the Air Force Geophysics Laboratory. 

The following is a description of a computer program developed for the CDC 6600 
computer which uses this technique. With the addition of a suitable user- 
supplied subroutine, the program can solve a wide variety of second-order 


ordinary and partial differential equations. 


METHOD OF SOLUTION 


A system of N coupled, second-order ordinary differential equations can 


be written in the following form: 


= 2 = a 
&' (x) = F(x) + A'(x) x F(x) + B'(x) B(x) = R(x) 
x 


Where A', 8', and €' are NxN matrices and Ff and R' are N-dimensional vectors. 
To solve this system of equations numerically, we first change to finite 


differences, Letting x =ndx (n=1,2,...): 


and 
| 
a2 ae net 7 2 Fo + Fy 
dx- 8x? 
the new finite difference equations become: . 
* R Fl 4 Bn Fh wi Ch Fuel . R, 


~ eK Rae) 
A = n iN n 
n 2 2 Ox 


| (5x) 
iy , SRO. 8 
: B= + B'(x) 
| n (éx)* n 
* C(x.) i A(x.) 
n (8x)* 2 Ox 
R, = R' (x) 


Boundary conditions at x) and Xy are assumed to be of the form: 


2) Oo) & Fl + Rv (x,) Fox) = RC) 
3) O(xy) LF, + AOy) Foy = ROD 


where x} is the value of x at the lower boundary and xy is the 
value at the upper boundary. In finite difference form, these 


are; 


we 


Ig 


nage oG) ~ . 2) 
Ry = 8 (a)) oe hy te 
C(x, C(x) 
8) =a ay Say + AD 
R = Ri (x,) Ry = R' (x) 


The solution (Richtmyer, 1957) is as follows. Assume: 


6) F =a Ff B 


+ 
n n n+l n 


where a, and BT are to be determined. Then: 


7) Fel = Ohl Fn - Biel 


Equation 7, when substituted into 1, yields: 


A, a By PR + Cf 


n n’ ‘n+l ~ (R, as A, B 


n-1 
Comparing to Equation 6 yields: 


f 

8) a =-(A a , +8 )-12 

| 

| 9) B= a, +87 mR, - A, By) 
At the lower boundary, using Equation 6: 


F, * a, F, + By 


Therefore: 


10) a, = -(A,) B, 
e 
11) B, = (A) R, 


Likewise, Equations 5 and 7 can be solved to give: 


= = ee = 
12) Fy = (Ay yey * By) Ry - Ay By) 
We now have the means of solving the equations. 


To review, we first calculate the finite difference form of the 
equations. a, and By are computed. Then we generate the other 0. and Bo 
using Equations 8 and 9, through a and By (note: oy is not needed and 
E, = By). The a and B, with Equation 6, generate the Fs completing the 


solution, 


PROGRAM 


The computer formulation is straightforward (see Appendix I for listing). 
The main program, TIDE, calls subroutines ABCR, ABCRN, ALPBET, and SOL which 


perform the steps outlined in the previous description. 


Subroutine ABCR is to be supplied by the user, This program calculates 
matrices A', B', and G' and vector R' for a given value of x. By proper 
manipulation of this subroutine, a variety of ordinary and partial differential 


equations can be solved. 


Subroutine ABCRN calculates the matrices K, B and C. and vector R, 
in the finite difference forms given in Equations 1, 4, and 5. Note that the 
matrices A, B, and C and vector R are reused. As the typical dimensions of 


such matrices are 40 x 40 or greater for our applications, this is a necessary 


process, 


Subroutine ALPBET calculates a, and BT from Equations 8, 9, 10, and 11. 


STREET 


They are stored on TAPE 4. TAPE 3 is necessary as a work tape for otherwise 
twice as many matrices would be required. If the matrix dimensions permit, 
this tape could be deleted and replaced by storage matrices, with substantial 


savings in computer time. 


Subroutine SOL uses Equation 7 to obtain the final solutions which 
appear in FNC. The output, in FNC, is printed, but it could be stored on a 


tape, plotted, or punched. 


EXAMPLE 


Appendices II and III contain the results of a sample calculation for 3 
ordinary, second-order differential equations. Appendix II is the exact 


solution while Appendix III contains the numerical solution. 


The example was defined as follows: 


0 0 ix ce 20 
A(x) = |o cael B'(x) = |0 0 0 
0 0 x? 0 oO x? 
2 0 oO (4-3 x°) + i(4 x°) 
Sruxjy =e fix® x2 0 Ri(x) = [its x2 + 3 x4) 
0 0 0 (4 x + x°) + (3 x? + 3 xi 
i=z/J-T_ 


The boundary conditions are, for x = 0 and x =T; 


* 
$ 
3 
a 


A'(0) = 


Oo 


on 
=> 
o 
— 
u 
Oo 


Pl 
oo 
J 
— 
" 
o 


= {0 


x4 + (3 x)i 


€*(T) 


B' (0) 


R' (0) 


B'(T) 


R' (T) 


3i 


u 
So 


aaa 


tT +3Ti 


i(6 T? + T) 


are: 
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C MAIN PROGR&M TIDE 
* a 
3 ~~ 
1 C PULPUS S 
t E NUMERICALLY INTEGRATES EQUATION CID2/0X2) F+eA(D/ OXF eOFar, 
c WHR ’e AgdyS AR® MATRIC S AND FoR ARE VESTOPS. AyBoCyP 
3 E AXE SIVEN 3Y SU3SROUTIN; ABCR. F IS TO BF COMPUTED 
' S 
0 UsAGE 
& PR)VIDI INPUT SUBROUTINE ARCR. OUTPUT AT STATTMENTESE 10 
Cc ANU 11. 
C 
Cc JFSCRIPTION OF PARAMCACTERG 


mF 


A > INPUT MATRIX 

3 - INPUT MATRIX 

5 - INPUT MATRIX 

R - [IVPUT VECTOR 

4¢T - STORAG: VECTOR 

FNUL - STORAGE VELTO2 

FNG = OUTPUT VESTOR 

BETecRM = DETERMINANT VALUE FROM MATRIX INVERSION. ILL- 
CONDITIONED PROBLEM WILL GIVE 0. 
IR = NUMBER OF 20WS 

IT - NUMBER OF SOLUMNS 

$1 - INITIAL VALUE OF X 

SN = FINAL VALUS OF X 

N - NUMBER OF STEPS 


REMARKS 
USER SHOULY CHANGE DIMSNSTONS AS KPPROF2TATE. IN THIS 
SAsE WHERE THE NUMBER 3 APPEARS. OUTPUT IS aT STATEMENTS 
19 AND 11. 


SUBROUTINGS AND FUNSTION sUBPROGRAMS FE QUIF ZY) 
AbCR - SETS UP PROBLEM 
APCRN = SONVERTS UsER?s PRIBLEM TO FINIT® DIFFERENCE FORM. 
ALPGET - COMPUTcCS INTESMEDIATE SOLUTIONS ALPHA ANO BETA 
SAIS - SYSTOM FUNCTIIN 
SOL = DOMPUTES FINAL SOLUTIONS F(X) 


METHUO 
SEC LINDZEN ANU KUD, MONTHLY WEATHOR REVIEW, VOLe S75 NO’ 
luy OST 1964,752-734. 


eeoeeveeveveeeeeoeveeveeeeee ev eeseeeeeeeereeeseeeoveaevneeoeeeseseaeseeaeee eee Beeeeeeeos 


OHO OMOOMOVGAAARH A OOA OKA Oe OOO L267 see 


FROGRAM TIDE CTNPUIL ,OUTPUT, TAP:i $2513, TAPES) 


C 
u CHANGE JIMENSTONS 
u 
COMPLEX AlSeg3) od C3 9509S 59309 05) 
COMPLEX 5° TO3) gFNE (3) 97 NO1L03) O56 TERM 
IR=5 
IC=s 
C ; 
% INTEGRATS FROME 


£1=0 


I u 103 
©N=2 
1; 
5) C IN STEPS OF3 
C 
Cx=.01 


NUMBER OF STESS WILL 328 


oO, 


R=C(SN“S1)/0X+1.001 
LOM FORMAT C2C1Ky -135.0)) 

bo 7 T=1,N 

¥=DX*#FLOAT(I-1)+#S1 


SET UP MATRICES 


30) 63 


CALL ABOR(A,3,CoR,TyNyXeIRyIS) 


COMPUTE MATRICES Ay3e9G AND VECTOR 28 IN FINITC ITFFEFPENCE FORM 


re eS) 


CALL ABERN(A,3 sColeNoI2R eo Io 9 1X) 


COMPUTE ALPHS AND SETA MATRT IS 


c> a 


CALL ALPSET(A,3yUy Re BHP yg ENClLs LyNo LT 2eICyDETERM) 
LET=Ue Bs (DETERM) 
IF(UFT) 7539397 
7 CONTINUE | 
COs T= ,2R | 
FNUC1)=B871 CD) 
& FNCL CI) =RET(T) t| 
19 FRINT 10G,FNO 4 
[0 5 J=2,N 
EACKSPACE 4& | 
BACKSPACE + 
FOAD(G) Ay,BET 


€> €3 82 


COMPUTE FUNCTLON 


CALL SOL (4,]5T sFNUgFNUL yIQy lV) 
11 FRINT 1ud,FNS 
fy LO o Url ik 
» FNCLCSSFNCCS) 
69 TO 97 
2 PRINT 100,09ET=RM 
‘7 CONTINUE 
STOP 
PND : 


RES AN NN RE OG EA ai eR a ce 


GOoOOOOOGCOPTS OOS Ste 


amos 


co 


G2 0a C2 C3503 62 07 O1E 


CHL es 


wn 


DULERJUTINE AZ3CR2 
PURPOST 
Usk a=P OVIIED GIGRYUILUF THAT SETS Ur FQUATICNS TO Ff 
SOLVLO. LOWER BLIUNDARY CONDITIONS ARE GIVEN IN STATEMENT 10 
TO Li, UPPER BOJNUARY CONDITIONS IN & TO 12, AND FOLUATIONS 
IN >» TO ve 


USAGE 
SET UP MATRICES 458,60 AND VECTOR R. 


DOESCLIPTION OF PARAMETERS 
AL - QJUTPUT MATRIX 
31 - OUTPUT MATRIX 
C1 = QUTPUT MATRIX 
Ri - OUTPUT VECTOR 


REMARKS 
USTR-PROVIOED. THE SOLUTION TO THIS EXAMPLE IS8 
FUKDSCX PFO, TE (KF FS) gh * G44 SFXFT) 


SUBROUTINES ANOQ FUNCTION SUBPRIGRAMS PEQUIFE) 
NONE 


METHOD 
USER SUPPLIED 


@esevneeeevevoeeevneeveen teevueeeeeeeaes eet eevee ovnevseeseev ees eee eeeeeee eevee eeeeesen 
SUBROUTINE AGORUAL »b1,C19R1xXsN,XoTR,IC) 
CHANGE UDIMENSTONS 


COMPLEX A1L(3,S)9B1L 035455) CL (55 5),R1 (3) 


1 I=1,18 
1 J=1,10 


AL1(I,J) =81(1, 3) =C1t (IL ,U) =CMPLX (0.059.0) 
LF (K-1) 259253 

C21 L 91) =01 029 3)YSCMPLX(Le ode) 
BL 2?y 3) =CMPLXU2053-59) 
£1°S$92)=SMPLK (329020? 

Ri lL d=RL(S)SSOMPLX(000,5.0) 
Fi1(2)=CMPLX(9.0,3.)) 

CO TI o 

TF (KAN) Soh4e% 
A1151)=CMPLX(105,0-0) 
L1¢3592) =UMPLX(V2Ze 95 eu) 

1 l2e3V SC MPLX (1069-9) 
h=X*xX 

b1(03592) =CMPLX(A,0.9) 

{=A*h 

Bl=567X 

CTHu oF K*X+A¥AtX 
fidi)=CMPLX(A,0.0) 
F1i(2)=CMPLX(B, BI) 
f£1(3)=CMPLXO0.0,CI) 

CO TO 6 


er 


| 


S C1191) =CMPLK (2005020) 
T1=X*xK 


Ti =T4*X 
C1251) S01 01,5) =CMPLK(0 20,72) 
C102 92)5A1 0252 FAL (3,5) =B1 03, 5)=0MPLX(T150.0) 
Fa4e-32*T2 
FI=4.*T5 
FLIC(1)=CMPLX(RYRID 
kI=c et T2+5 ete 
1 (2)5CMPLX(9.0,FI) 
F=4.*T4+TS 
FI=5e*T1+3.*T? 
K1i(3)=CMPLACR,RI) 

& CONTINUE 
RETURN 
ENC 


8 SUBROUTIN= ABC2N 
is : PURPOSE 
: : CONVE21S OLFFE2ZNTLAL FORM OF EQUATIONS [9 FINITE CIFFER- 
| 2 ENCS FORM. 
& ~ 
j y 
c USAGE 
5 INPJT MATRICES AyByC AND X VALUE. OUTPUT IN AyBeCe 
v VESCRIPTION OF PARAMET=RG 
5 A - INPUT,OUTPUT MaTRIy 
5 3 = INPUT,OUTPUT MATRIX 


€ = STOP NUMSTFE 


u No - FINAL STEP NUMBE2 
C If - NUM3ER OF ROWS 
u Iu - NUMBER OF COLUMNS 
C OX = STEP SIZE 
e 
o REMARKS 
o DIMENSTONS MUST eF CHANG ID 
C 
o SUNYROUTINES AND FUNCTION SUBPROSRAMS REQUIKED 
Cc NONE 
C 
C MF IHOS 
C FINIT= OIFFERENSS APPROXIMATION TO DIFFERENTIAL. 
u 
Cio 6 O66 6 6.6.6 6 6 OO OHNO © 010.68 66 66.0066 00 00 60 600 40.5 04 6 OOOO OKC MOR CBE OSOC ROO 2 08 
C 
SUBROUTINE ABLCENCA gD els KoNeIRo IC 9 OX) 
C 
c CHANGE JIMENSTONS 
CG 
COMPLEX A393) 550555),50353) 
Fri./UX 
1F(K=-1) 151,52 
2 IF(K=N) S58,8 
1 (oO tu IT=1,1R 
CO 13 J=1.Id 
A(I, JV =ACI oJ) -F*FC (TL, JD 
10 E¢(Iyu) =F*FC(I sy) 
Co TO 6 
3 t 2=F +F 
F3=F/e6 


Fus=-ce*F2 
CO 2c I=1,1R 
e tO 2C€ J=1,Ic 
EdLs JI HBC J) + F4*C (1 _J) 
CIL,vI=F2*C(l, JD+F 3* ACT, eo) 
20 R(T JI=HCCI JI -FFACTI, J) 


CO fO 6 
3 CN 3uv 1=1,1R 
C0 $G J=1,1C 
UII, JI=ACT J +F*YC (ly J) 
30 ACL,» uv) =-F*O(T, Ud 
t RETUSH 
| 
' 
| 


ENO 


SUBROUTINE ALPOET 


nha 
rw) 


ee 


Rsa0e 


PUP POSE 
CaLCULATES ALPHA ANU 3¢TA WHICH Akt NEFNEL AS INTERMEDIATE 
SOLUTIONS. OUTPUT A®PPTARS ON TAPLG. 


USAGE 
PROVING INPUT MATRIC-35 Ag3ec AND VCCTOF 1. it TURN ALPHA 
AND BETA ON TAPE4G. 


OFSCRIPTION OF PARAMETERS 
A - INPUT MATRIX 
3 - INPUT MATRIX 
CG - INPUT MATRIX 
0 - INPUT VECTOR 
3ET - WORKING VECTOR 
FNC1 - WORKING VECTOR 
K = STEP NUMBER 


2OOQ OWE AGIQ|aAaIGAVYGAAeA 


~ 
‘ 


1 ae Pa Ob Oe BN eRe 


5 N - FINAL STEP NUMBER 
' c IF - NJMBFR OF 2OWS 
| 
c IC - NUMBER COLUMNS 
¢ DETIRM - VALUD OF DETERMINANT 
C 
3 REMARKS | 
C ThPE3 IS A WORK TAPE. IW MANY CASES IT CAN BE REPLACED ! 
C WITH MATRICES TO CUT RUN TIME 
3 SURROUTINES AND FUNSTION SUBPROGRAMS REQUIRED ) 
0 MINY = MATRIX INVERSTOU | 
S MPR00 - MATRIX MULTIPLICATION 
C GMaDD - MATRIX ADDITION 
C 
G MEY HOD 
L SEC REF ERENCE 
c 
Cb WE CE 66-6 OEE. 6 8E SEEK CEES ECOEECERE CEH EEM CECE STEN CHOOSES COCKS 6666 0 4180S OO 6S 
u 
SUBROUTINE ALPBET (Ay 34090, BIT y ENOL eKyNoIRy ICyDETERM) ) 
; CHANG? UITMENSTONS ‘ 
COMPLEX A(3—3) 980393) 95 (395) 5 003) BET (C3) 9F NGL (3) »DETERM ; 
FORMAT (6 (1X9=1306)) 


191 FORMAT (1X 51+) 
IF (Kei) 19is2 
1 CALL MINV(4,I° ,U51,0ETIO RM) ‘ 
(9 11 KI=1,12 
10 11 JI=1,15 
11 A(KIsuI) =-AC(KI,JI) 
CALL MPRLD(Ag' gCyIRy IU, IC) 


— 
a 
=) 


Pi Sart te siti 


2ITi ALPHA ZERO AND BATA Ze 20 ON TAPES S ANN G 
TaPFo IS A WORKING Tare, TAPS. STIRTS ALL ALPHA N?S ANG BETA NAS 


| fe . 


14 


VOITECS) Usd 
WRITE (G4) Sy,0 
F FETURT 


CONTINUZ 


2 
} 
c fTORe CN AND GN ON T4ris 
U 
bRITE(3) Cy3 
FEWIND § 
C 
G FEAU WLAIN=1) ANDO SETCN*1) F2AOM TAPT3 
i G 


READCZ) Cy SET 


CALL MPRDD(O,8ET,FNCLeIr 1Cyl) 
CALL GMACOCU,FNC1I,BETsIr oly -1.)? 
CALL MPRUD(M,0,939,T Re IU,1C) 
CO & IT=1,IPR 
CO 4 J=1,IC 

4 ACI, J) =H8(T J) 


FeTXI:VE ON AND BN FROM TAPES 


FFAD(S) 0,3 
CO S LT=1,Ik 
{QO 5 Jsi,Iec 

5 C(I,y=-ClIyJ) 
FEWINU 3 


CALL GMAHD(Bs4,3eIReIU, 16) 
CalLL MINWVCoeTl p3ET aly IE TERM) 
CALL MPRNU (Bs, .,4,1R,IC, 10) 


3 VRITe BETOND GNI ALON) ON TAPFS 
C 
heIde (3) ABET 
C 
c WRITE BCT(N) &NO ALON) ON TAPE Ge 


wEITE CG) Ay3iT 
FETUPN 
Fo 


om 


Co 


SK oy <p Mor a op ew We OM Bw Gk Or ee op Je 


ae 


CF C253: 47 82.07 62 


55 


Pe! 
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SULROUTEINE SOL 


PUPPISe 
GIVEN ITNTERMEDJTATE SOLUTIONS ALPHA AND 8£TA AND PREVIUUS 
VALUE OF F(X), SOLVES FOR FUX-UX). 


USAGE 
PROVIOF ALPHA, BETA, AND F(K). RETURNS FUX-2PX), 


DESCRIPTION OF PARAMETERS 
ALP - INPUT MATRIX ALPHA 
3ET - INPUT VeEuTOR BETA 
FNCL - PREVIUUS VE STIR SOLUTION F(X) 
FNG - XETURNED VECTOR OLUTION F(X-0X) 


ROMARKS 
NONE 

SUPRIUTINES ANU FUNSTION SJUBPROIOGRAMS REQUIP=ED 
GM4AD9 - MATRIX ADCITIIN 
4PRIO - MATRIX MULTIPLICATION 


MeTHON 
SFE REFERENCE. 


@eeeceeoevoeneesneeveeveoeee eeeoeeeeeeeveeeseeoeeeeceoeeteeseeevneseseeeeeeeeeeeeeoeeeee ees 
SUBROUTINE SOL CAL?,3EF,FNOyFNU1L,IR,IC) 


CHANGE JDJIMENSTONS 


COMPLE X ALPO593) 5 3°5T (35) gy FNO(3) gFNCL(S) 
CalLlL MPAIVD(ALPyFNUL FNC sIR,ICg1) 

CALL GMASUCFNG »BET FNC, Ikyi,yi.?) 

RETURN 

END 


SUBROUTINE MINY 


PURPUSF 
MATRIX INV: RSTION WITH ACCOMPANING SUOLUTTONS GF LINEAR 
cQUATIONS. 


USAG 
INPUT MATRICES A AND 36 RETURNS AS INVOFS IN A AND ASINVERSE 
TIMcS §} IN 8. 


DESCAIPTION OF PARAMETERS 
A - INPUT MATRIX TO 3f INVERTED 
INVERS® RETIRNED Til Ae 
3 - MATRIX OR VECTOR 3UCH THAT AX*5=B. 


Cr) 0929. Ca £7 09 Ce Oo 6O Ce 


we 


C A-INVERSE TIMES 2 FETURNED IN 3. 

o N - SOLUMNS IN Ay ROWS IN 3. 

v M - SOLUMNS IN 3. 

C JETER = DETERMINANT OF A 

Cc REMARKS 

G Fi N&S3 MATINV MATRIX INVERSION 

o SUBROUTINES AND FUNSTIION SUBPROGRAMS R© QUIRED 

c CABS = SYSTEM FINCTIOIN 

C 

o METHOD 

C GAUSSTAN ELIMINATION WITH PARTIAL PIVOTING. 
Lv 

Cre ee ee ee 0s rr rs 
c 


SUBROUTINE MINV(A, Ny 394, 0ET=2 4) 


C CHANEL DIMENSIONS 
{ S 
CIMERSTION 806353) 55(3 51) ,PIVIOT(3) »IPIVOT(3) » INDEX (542) 
COMPrLS x Ng 3 JETERM,PTVOT sAMAX ySWAP gT 
C INITLACI2ZA7TIIN 


19 LETExh= CMPLX (1.2050.9) 
15 DO gu J=1_N 

20 1PIVvul (J) =0 

ju 00 Sot T=1,_N 


L 
Cc SCARCH FOR PIVOT <LEMENT 
C 


. 4&9 BMAX= CMPLX (G+Ugued) 

45 £0 105 J=49N 

50 IF CIPIVOT(J)-1)00,195,60 

bO £€O 15u KELyN 

70 IF CIP IVOV(K)-19305 1099740 

30 TF CC AL SCPMAX) CABS (A (dy 899) 3591905100 / 
85 = 1°0N=u ; 
ee) TLOLUM=K j 
35 6 AAX=B Sy) 
100 CONTINUS 

105 CONTINUE 


17 


eareren maaan: 
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AJL TWO GENt RAL MATRIC-% TO FORM RESULTANT GENEFAL MATRIX 
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CALL bM6 0D 04 43,2 aNoM 9A05) 


ASS CR OF PARAMETERS 
~- NOME OF FIRST INPUT MATRIX 
s ~ NAM® OF SESONU INPUT HATRIX 
R = NAMES OF OUTPIT MATRIX 
N ~ NUMSER OF ROWS IN 44852 
M - NUMBER UF COLUMNS IN 4,3, 


REMARKS 
ALL MATRICES MUST BE STORED AS GENFRAL MATPTCES 
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SUBRIUTING MPROD 


PURPOSE . 
MULTIPLY TWO GEW2RAL A°VIRIVES TO FORM A 2PFSULIANT GENERAL 
MATRIX : 

USAGE 


Call MPOD CAS 9 co Noel) 


CESCRIPTLON OF PARAMETERS 
A - NAM” OF FIRST INPUT MATRIX 
B - NAW OF SESCONU INPUT MATRIX 
R - NAM® OF OUTPJIT MATRIX 
N - NUYSFER OF RUNS IN A 
M - NUYBER OF SOLUMNS IN 4 AND ROWS IN & 
L - NUYRER OF COLUMNS IN 8 


REMARKS 
ALL MATSICFS MUST BE STORED AS GENERAL MATRICES 
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MATIX * CANNOT 3€ IN THE SAME LOCATION AS MAITPIX 8 
NUMBER CF COLUMNS OF MATRIX A MUST 36 FLUAL TQ NUMBER OF ROW 
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SUBROUTINES AND FUNCTION SUBPROGRAMS PRIQUIREN 
NUNE 
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THE M BY L MATRIX & IS PREMULTIPLIED OY THE N BY M MATRIX a 
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